A KKR Green function formalism for ballistic transport 
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We develop a method for the calculation of ballistic transport from first principles. The multiple 
scattering screened Korringa-Kohn-Rostoker (KKR) method is combined with a Green's function 
formulation of the Landauer approach for the ballistic transport. We obtain an efficient O(N) 
algorithm for the calculation of ballistic conductance through a scattering region connected to semi- 
infinite crystalline leads. In particular we generalize the results of Baranger and Stone in the case of 
Bloch wave boundary conditions and we discuss relevant properties of the S-matrix. We consider the 
implications on the application of the formalism in conjunction with a cellular multiple scattering 
description of the electronic structure, and demonstrate the convergence properties concerning the 
angular momentum expansions. 

PACS numbers: 72.10.Bg,72.15.-v,71.15.-m 



I. INTRODUCTION 

The study of transport properties of solids has always 
been a challenge for experimental and theoretical con- 
densed matter physics. The difficulty lies among other 
things in the fact that the transport properties are by 
definition related to non-equilibrium situations, so that 
their modeling is not easily based on standard techniques. 
Fortunately, in the regime of linear response to weak ex- 
ternal fields, the study can be based on ground state 
properties, by treating the field as a perturbation. That 
this approach is founded solidly is guaranteed by the 
fluctuation-dissipation theorem; in particular for the case 
of electrical conductivity this is expressed by the famous 
result of Kuboi 

In the past decades, the development of first-principles 
methods for the calculation of the electronic structure 
of solids has been accompanied by corresponding ad- 
vances in the formalism and methods for the calcula- 
tion of transport properties. Technological interest has 
given a push to the field, and novel effects such as the 
Giant and the Tunneling Magnetoresistance (GMR and 
TMR) have been investigated theoretically and applied 
in technology. Furthermore, the size of today's elec- 
tronic devices is so small that the understanding of bal- 
listic (i.e. phase-coherent) transport has become impor- 
tant not only for basic physics but also for applica- 
tions. As a result of these developments, ideas and meth- 
ods which initially were conceived for the understand- 
ing of electronic transport in simple cases, have been 
combined with techniques based on a realistic descrip- 
tion of the electronic structure in order to give reliable 
and material-specific results. 2 Without claiming to give 
a complete list, we mention that for diffusive transport 
through disordered systems there exist methods and re- 
sults based on the combination of the Korringa, Kohn 
and Rostoker (KKR) Green function method with the 
Boltzmann formalism^^ on the coherent potential ap- 
proximation (CPA) combined with Kubo-GreenwoocU^ 
theoryjiSi 9 . or on the tight-binding method^ For bal- 



listic transport there exist methods based, for example, 
on LEED; 11 ! 12 layer-KKRi^. or similar layer-type^ tech- 
niques, on the tight-binding approachfi^iS*ii and on the 
transfer matrix concept ^I&ia these mostly combine the 
Landauer-Biittiker approachS2i2i with electronic struc- 
ture methods, as will be done in the present paper. More- 
over, due to potential applications in GMR and TMR 
devices, spin-dependent transport has come to the cen- 
ter of interest with emphasis to conduction in magnetic 
multilayers 2 *£*2i±&i22i22i2i and ferromagnet-semiconductor 
hybridesi 1 ?' 1 ?: 2 ^ 2 ^ 2 ?! 2 ?:??!? 1 In these systems the elec- 
tronic spin degrees of freedom are accounted for in or- 
der to achieve spin-dependent resistance. In addition, 
novel systems such as nanowires or atomic-size contacts 
are created experimentally, and demand interpretation of 
their transport properties. 

In this paper we present a method for the calculation 
of ballistic transport from first principles, which com- 
bines the KKR ab initio Green function technique with 
the Baranger and Stoned formulation of the ballistic 
transport. The method is suitable for layered systems 
and interfaces with two-dimensional periodicity, as well 
as for atomic size constrictions connecting two leads or 
nanowires^ and supports spin dependent effects. Since 
the KKR technique in the screened formalism offers lin- 
ear scaling of the computational effort with the number 
of layers^ our method can handle large systems. Ad- 
ditionally we present the proof of some theorems con- 
cerning the transmission probability of wavepackets in 
crystalline, rather than free-electron, environment. 

The paper is organized as follows. In Section[n]we give 
a description of the setup of the physical systems that our 
method can be used for. In Section ITU we briefly address 
the approximations and assumptions made, also in con- 
nection to the Kubo-Greenwood and the Baranger-Stonc 
formalism. Sections IIVI and [V] are devoted to making 
the connection to the Landauer-Biittiker formalism, and 
some elements of scattering theory are given there. The 
conductance formalism for KKR is developed in Section 
IVII Some examples illustrating the convergence proper- 
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FIG. 1: Setup of the problem. The regions L and R corre- 
spond to perfect semi-infinite crystalline leads (not necessarily 
of the same material), attached to an "interaction region" I. 
This includes the interface or other structures sandwiched be- 
tween L and R, plus a few monolayers of the leads at each side 
so that the evanescent interface states can decay. The con- 
ductance is evaluated between the surfaces Sl and Sr., where 
only bulk Bloch states should exist. 



ties of the method are given in Section fVIII Finally, we 
conclude in Section IVIIII The appendices contain parts 
where lengthy mathematical manipulations were needed. 



II. SETUP OF THE PROBLEM 

The systems that we study consist of two half infinite 
perfect crystalline leads, left (L) and right (R), attached 
to a slab which is considered as the "interaction" region 
(I); schematically this is represented in Fig. ^ The sur- 
faces iSl and «Sr should be chosen far enough within the 
leads so that any localized interface states have decayed 
and only Bloch states are present there. The direction 
of growth is taken to be the z direction. We consider 
systems with a 2D periodicity in the xy plane where the 
scattering region is embedded between two semi-infinite 
crystalline leads. Moreover, summing up current contri- 
butions in real space we can calculate the conductance in 
the presence of defects in nanowires^i^ and we can even 
treat atomic sized constrictions between infinite leads if 
the current flow is localized in the constriction region 
which is usually the case44 Thus we can simulate trans- 
port through small molecules or a break junction geom- 
etry. 



III. KUBO FORMALISM 

The connection between the ballistic dc conductance 
or conductivity and the one-electron Green function 
in the linear-response regime has been given in the 



past 



32.35,36.37 



In short, one uses first-order perturbation 



theory to calculate the effect of a weak, time-oscillating 
electric field on the density matrix of the system; then, 
the frequency of the oscillating field is taken to zero, and 
the dissipative term of the expectation value of the cur- 
rent gives the dc linear response of the system. The re- 
sult is just an expression for the Kubo conductivity of 
the system at the limit of zero temperature. Integration 
by parts of the conductivity expression gives the relation 



between current, conductance and voltage. Baranger and 
Stoned have proven that the conductance is a Fermi- level 
property, although the two-point conductivity within the 
sample can also have contributions from other energy lev- 
els when magnetic fields are present. Here we briefly de- 
scribe the relevant formalism, in connection to the rest 
of the paper. 

The general relation between the current density at 
some point, j(r), and the electric field at all other points 
of a material, E(r'), in first order perturbation theory for 
a stationary state, is given by 



JW-/ 



= / dVer(r,r')E(r') 



(1) 



The nonlocal conductivity tensor er(r, r') is connected to 
the retarded one-electron Green function G + {v,v' \E) of 
the system via the famous Kubo-Greenwood result. Its 
frequently used form is 7 



<x(r,r') = --J dE(-f'(E)) (2) 

x ( —VSG + (r,r';E)) ( -V'SG + (r',r;£)) 
\m J \m J 

Here, f'(E) is the Fermi function derivative and 3 de- 
notes the imaginary part. In terms of the difference be- 
tween the retarded and advanced Green function G~ 

AG(r,r';E) := G+(r, r'; E) - G"(r, r'; E) (3) 
= -2m [ dai> a (r)i>* a (r')5(E - E a ) (A) 



the conductivity can be written 



ff(r.r') = - 



1 ft 3 



dE(-f'(E)) 



16irm 2 

xAG(r,r';£)VV'AG(r',r;E) (5) 

< — > 4 — > 

where the symbol V is defined as /(r) Vj(r) = 
/(r)Vg(r) — (V/ (r)) g(r). Here only the delta- function 
part of the Green function has been kept; the principal- 
value part gives rise to the reactive term, which is zero 
for dc conductance^ Upon taking the limit of zero tem- 
perature, the derivative of the Fermi-Dirac function be- 
comes a delta function according to —f'{E) — > 8{E — E F ) 
and contributions to the conductivity only arise from the 
Fermi level. 

In the absence of a magnetic field, as assumed in our 
work, the nonlocal conductivity tensor is symmetric un- 
der the interchange of r and r', and divergenceless: 

<x(r,r') = <x(r',r) (6) 
V-er(r,r') = <r(r, r') • V' = 0, (7) 

where V means that the operator acts to the function 
on its left. The former relation is a consequence of the 
symmetry G(r,r') = G(r',r) of the Green function for 



3 



local potentials. The latter is a result of current conser- 
vation in the absence of magnetic field; if a magnetic field 

is present, one has merely V • <r(r, r') • V = instead^ 
In Eq. (JIJ, of course, E(r') is the total field, i.e. the 
external one plus the one being induced by charge relax- 
ations; the electronic gas and nuclei must be also fully 
taken into account. This is not easy to handle; however, 
this expression can be integrated over the whole sample 
leaving only the current in relation to the applied voltage 
by exploiting the fact that E(r) = W(r), where V^r) is 
the electrostatic potential. Inserting this into Eq. Q one 
gets for the total current J through Sr 

J := f dSzj(r) = AV I dsf dS'za{r, r')z', (8) 

after an integration by parts and use of Eq. J7J). Here AV 
is the external bias voltage. The "zero-current theorem" 
stating that the ground-state electrostatic potential gives 
always zero current density^ has been used, so only the 
externally applied voltage AV remains. Thus one can 
recognize the conductance as the flux of the conductivity 
tensor through the surfaces Sl and Sr: 



9 



[ dS [ dS'z<r{v,v')z' 

J St, J Sr 



(9) 



Substituting Eq. JSJ) one obtains for the conductivity: 



2^3 



fJ 



dS / dS' 



167rm 2 7 Sl 7 5r 
:G+(r, r'; E F )V 2 V ' 2 G~(r', r; E F ) (10) 



where the terms involving G + V z V ' Z G + and 
G"V Z V 'G~ vanish^ 



IV. LANDAUER FORMALISM 

In the Landauer approach, the conductance problem is 
viewed from the aspect of scattering theory. In this way, 
the causal relation between voltage and current is con- 
ceptually reversed^ Instead of applying a voltage and 
examining the current as the response, a current is forced 
to flow through the sample and the voltage is viewed as 
the result of the pileup of carriers at the various obsta- 
cles, forming residual resistivity dipoles. The result is, 
of course, again the usual current-voltage relation. In a 
multi-lead experiment, if each lead n is held in potential 
V n , then the current J n going out through this lead is 



(11) 



The coefficients g nm describe the conductance of the sys- 
tem, and in the case of only two leads there is only one 
conductance coefficient g. 

The concept of incoming and outgoing scattering chan- 
nels is introduced, which play the same role as in- and 



out-states in scattering theory; in the cases of our interest 
they are Bloch states in the leads. We need to describe 
the scattering process of one-electron Bloch states inci- 
dent from L as incoming waves and scattered into L or R 
as outgoing waves; thus an S-matrix formulation is ap- 
propriate. One has waveguide or crystalline geometry in 
the leads, rather than the free-space geometry of usual 
scattering theory, and in addition the leads can consist 
of different materials. In this respect, usual scattering 
theory needs a few modifications to be applicable. 

Once the S'-matrix elements are available, one can 
readily calculate the transmission probability T„j ;rra j 
from each incoming channel i in lead m to each outgoing 
/ in lead n, and use the Landauer-Buttiker formula to 
calculate the conductance g nm from lead m to lead niSi 



/ nf;mi 
fi 



= tl> 



nf;mi \ 



(12) 



fi 



Stone and Szafer— and Baranger and Stoned showed 
the connection between the Landaucr-Biittiker formula 
with the Kubo result starting from perturbation theory, 
considering free-electron leads. In Sec. we shall follow 
their analysis closely, but prove some of the theorems 
needed for crystalline leads taking into account the Bloch 
character of the incoming and outgoing channels. 

The S'-matrix formulation in such problems and the 
question of unitarity will be shortly addressed now. First 
of all, the in-states and out-states, which are in a usual 
formulation plane or spherical waves, are here propagat- 
ing Bloch states in the leads. Whether such a state 
is incoming or outgoing is determined not by the Bloch 
wavevector k, but by the group velocity Vk = Vk-Ek- 
Suppose that a particular lead is grown in the z-direction, 
with the unit vector z pointing away from the sample. 
Then, the Bloch state is outgoing if (vk) z > 0, and in- 
coming if (vk) 2 < 0. Evidently the set of in-states can be 
different than the set of out-states. Consider, for exam- 
ple, the case of only two leads of different materials, say 
L in the left with wavefunctions 5 ,L and R in the right 
with wavefunctions , J R , sandwiching the interaction re- 
gion I. The set of in-states for this problem will consist 
of right-traveling waves and left-traveling fy^, while 
the set of out-states will consist of left-traveling waves 
^o Ut and right-traveling waves ^f^ ut ■ 

The total wavefunction in the system will be asymp- 
totically a linear combination of in- and out-states; in 
particular, if we choose incidence from L as a boundary 
condition, we have the form 



*k» = 



ka 



Ek'a' Hcak'a'l'k'a"' 



E, >TrRout 
k'a' t kak'a'y'k'a' 



-00 



(13) 

For finite z values also evanescent states have to be 
included in the summations, however they die out for 
z — > ±oo. Here, band indices a and a' and Bloch 
wavevectors k and k', have been introduced; elastic scat- 
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tering is implied. The total wavefunction (r) is char- 
acterized by ka in the sense of the boundary condition, 
i.e., it originates from an in-state which has this wavevec- 
tor. The transmission amplitude ikak'a' and the reflec- 
tion amplitude ncak'a' are related the elements of the 
S-matrix of the system. 

The normalization of the in- and out-states determines 
whether these amplitudes t and r are identical with the 
S- matrix elements or not, since the S- matrix must al- 
ways be unitary. If one normalizes the in- and out- 
states to unit flux (as \P — » ^/y/v), then the transmis- 
sion probability from state \t£ a in to state V^'a' out is just 
Tkk'a';Lka = |<kak'a'| 2 = IS'kak'a'l 2 - But with the usual 
normalization to unit probability in space, the transmis- 
sion probability must invoke the group velocities to ac- 
count for the different flux of in- and out-states^ 

Tkk'a'jLka = |<5kak'a'| 2 = |tkak'a'| 2 77 (^) 

\\ v ka)*\ 

This point will be addressed in detail in the following 
section. 



(^ka ano - ^k"!' nere are supposed to be left- and right- 
traveling states in the same lead). These have been 
proven before^ for free electrons; we present the proof 
for Bloch states in Appendix 1X1 

Eqs. (|15I16|> can be used to project a scattering wave- 
function onto a particular channel; using them in con- 
nection with Eq. 1)13(1 we can extract the transmission 
amplitude out of 



9m C <— > 

(17) 

(it) Next we consider the asymptotic expression for the 
retarded Green function G + (r, r';E), with r and r' in dif- 
ferent leads and going to infinity. This can be found^S 
by using the boundary condition that G + (r, r'; E) should 
represent an outgoing wave at r, G~(r, r'; E) should rep- 
resent an incoming wave at r, and 



V. CONNECTION BETWEEN KUBO AND 
LANDAUER APPROACHES 

The connection of the Landauer approach involving 
the S-matrix to the Kubo-Greenwood conductivity for- 
mula involving the Green functions has been given by 
Fisher and Lee2&, Stone and Szafen^i, and Baranger and 
Stoned in the case of free electrons in the leads. How- 
ever, for crystalline leads one must account for the rele- 
vant band structure. Here we present some proofs needed 
to extend the above results to crystalline leads. In par- 
ticular we shall pursue the S-matrix elements between 
in-states from L and out-states in R and see how these 
connect to the Green function of the system and to the 
conductance. In this way the connection of the Kubo 
to the Landauer formalism will be made. We proceed in 
three steps: (i) we find an expression for the transmission 
amplitude ika;k'a'! (H) we express the asymptotic Green 
function in terms of the S'-matrix; and (in) we express 
the S'-matrix in terms of £ka ; k'a'- 

(i) We start with the calculation of current matrix 
elements between Bloch states. In particular, let ^g 1 
and ^^a* be incoming and outgoing Bloch states (in the 
sense described before, i.e. right- and left-traveling) of 
the same energy (E^ a = Ek.'a') in the same lead; the 
states are normalized to unit probability rather than unit 
flux. Also, let S be the lead cross section, normal to the 
z-direction. Then the following orthogonality relations 
hold for the current matrix elements: 

/ ds*k* vUi?w = / dsnT v»*e& 

Js Js 

= ■i—\(vka)z\Skk'Saa' (15) 

/ dS*&* V Z *S = (16) 
Js 

(for £ kv = E ka ) 



G + (t,t';E)=G-*(t',t;E). (18) 



The Green function then, expanded in in- and out-states, 
has the asymptotic form 



G + (r,r';E)^Y,T, A ^'^a Ut (r)^r(r') (19) 

ka k'a' 

(for z — > oo, z' — > — oo) 

We are now seeking a relation between the coeffi- 
cients ^4ka;k'a' and the transmission amplitudes £ka ; k'a' ■ 
Eq. I|17|) was used together with the Lippmann-Schwingcr 
equation and Eq. (f 1 Of) in Rcf. 32 to prove this for 
free-electron leads. However, the Lippmann-Schwingcr 
equation is difficult to handle when we have differ- 
ent materials in the leads, because it connects the 
Bloch wavefunctions of different materials in the infi- 
nite leads, thus there is no localized perturbation. Al- 
ternatively, one can start directly from the definition 41 
of the S-matrix in terms of the time-dependent retarded 
Green function (propagator), looking at the transmis- 
sion amplitude from an initial wavepacket $^ in (r, t) — 
J dEie-^a^Ei)^™^; E,), incident from the left, to a 
final one $^ out (r,t) = / dEfe'^afiEf^f^^Ef), 
outgoing to the right. Here ^f?i and if? f denote itinerant 
Bloch states and we have propagated the initial and final 
states to t and t' by their corresponding bulk hamiltoni- 
ans giving the exponential factors. At the end aj and cif 
will be taken as extremely peaked distributions around 
the same energy Ei\ i and / represent then definite k'a' 
and ka. One gets 



5 



Sfi := 



lim 
lim 
lim 



/" d 3 r / d 3 r'$^ out *(r,t)G + (r,i;r',t')^ in (r , ,i') 

/" d 3 r y dV$* out *(r,i) J dEe~ iB ^-^G + (r, r'; S)$}' in (r / , *') 

y d J B / e iE / i a}(E / )(^ out (E / ),^ out ( J B)) y d^e^'a,^^^),*^^)) 



= y dEa* f (E)ai{E)viVfAfi(E) 
= A f i(Ei)y/ViVf (on-shell). 



(20) 
(21) 



The velocities appearing represent the z component of 
the group velocity, since in waveguide geometry the k|| 's 
form in reality a very dense discreet set so that the 
wavepacket is constructed for a definite k|| from the con- 
tinuous k z spectrum. Having this in mind, we have used 
in the derivation of (|2U|I the orthonormality condition for 
Bloch waves 

(* ka (E f ),*v a ,(E)) := Jd 3 rn a (r,Ef)*Va'(r;E) 

= S^' Saa'Sikg - k' x ) 

= 5^ k ,6 aa ,v z 5(E - E f ). (22) 

In addition we have used the normalization of wavepack- 
ets to unit probability 



1 



dEa*(E) y dE'a(E')(^(E),^(E')) 
dE\a{E)\ 2 v E (23) 



which for a very peaked distribution a{E) around some 
energy E$ gives 



dE\a(E)\ 2 = l/v Ei 



(24) 



implying that in the limiting case |a(E)| — > 
y/ 8{E — Eij/vEi- Finally, in the last step we 
have assumed that aj{E) and a>i(E) are both 
peaked around the same energy Ei so that the 
wavepacket goes to a single Bloch function, whence 
a}(E)a t (E) -> 5(E - E,)/ ^v f {E l )v l {E l ). Thus, the 
final result 1)21(1 is valid for on-energy-shell scattering of 
Bloch waves; else, for general wavepackets, Eq. (|20f) must 
be applied. Working with wavepackets has guaranteed 
the correct normalization. 

From Eq. (|21|l we see that the coefficients in the Green 
function asymptotic expansion are just 5-matrix ele- 
ments (normalized to the group velocities). In this form, 



the S- matrix is unitary, i. e. the scattering probability is 



T 



fi 



\S 



ft\ 



(25) 



(Hi) Finally we show that the relation of the S- 
matrix to the previously defined transmission amplitude 
of Eq. JT3J) is 



(26) 



i.e., a normalization involving the group velocities is 
needed. This can be seen by noting that an incom- 
ing wavepacket $i n = J dEa(E)^/ m (E), normalized to 
unit probability as in Eq. (J2HJ, evolves partly into the 
wavepacket $ out = J dEt(E)a(E)^ out (E) according to 
Eq. I|13|) . Assuming that a(E) is so much peaked 
around Eq that t and v are constant in this energy 
range, the scattering probability is given by the nor- 
malization factor of the outgoing wavepacket: T/j = 
||$out|| 2 = \t\ 2 J dE J dE'\a(E)\ 2 ^ out (E)^ out (E')) = 
\t\ 2 J dE\a(E)\ 2 J dE'5(E - E')v out = \t\ 2 v oat /v in , where 
the normalizations 1221) and (|24|) have been utilized. 
Comparing this to Eq. proves Eq. (|26l up to a phase 
factor, which can be seen to be just unity by creating a 
wavepacket out of ^ tot in Eq. (|13fl and constructing the 
S- matrix element. Furthermore the expression i|14|) re- 
sults from (t2"5|) and 

Combining all the above, we may rewrite Eq. (|20|) as: 



G+(r,r';£) 



Ska; 



k'a' 



EE 

ka k , a , \/(vk'a')z(vka)z 



*^(r')*£, out (r) 



EE^f^(r')C'W (27, 



ka k'a 



7^ (Vk'fl')i 

' n' v 



which gives us for the advanced Green function 



G-(v,v';E) 



EE 

ka k'a' 



ka;k'a' ^Rout* lT( L 



(Vk'a' 



k'a' 



(r')*^ n (r). (28) 
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We can extract the ^-matrix elements from the Green 
function: 

h 2 1 

Ska;k'a> = ~. J . (29) 

xf dS< f dS^,(r')W' z ,G + (r,r';E)W z ^ a ° Ut *(r) 

where Sl and Sr have to be sufficiently far apart so that 
the asymptotic formula i|2(J|) without evanescent states 
can be used (the evanescent states are always included in 
the self-consistent Green function). 

As a final result, we combine the above steps (i)-(iii) 
to arrive at the Landauer formula. We start from the 
Baranger-Stone expression (fTUl) , expand the Green func- 
tions according to 127|) and (J2SJ, and use the orthogonal- 
ity relations (|15|l to get rid of some terms: 

g = -- j / dS / dS' 

16irm z J Sl J Sr 

xG+(r, r'; E F )V Z V ' z G"(r', r; E F ) 

ka m k'a' out 

= t S S Tk ' a ' ;ka (30) 

ka in k'a' out 



which includes the local t-matrix of each single-sited po- 
tential. For further details we refer to Refs. 133. 

The layered systems that shall be considered with the 
KKR formalism consist of two half-infinite crystalline 
leads, assumed to have perfect periodicity otherwise. 
Sandwiched between these leads is an "interaction" re- 
gion where a different material can be placed and where 
the scattering of the Bloch waves takes place. Three cases 
can be handled in this respect: (i) systems with in-plane 
periodicity^ (ii) wire-like structures embedded in vac- 
uum (or in some other non-conducting medium) and 
(Hi) atomic constrictions between semi-infinite leads. 

Concerning case (i), when we consider systems with 
two-dimensional in-plane (x-y) periodicity (perpendicu- 
lar to the direction of growth z), the interaction region 
and the two leads have common in-plane Bravais vec- 
tors. If needed, larger (non-primitive) two-dimensional 
unit cells are taken to match the lattice constants of both 
materials; this is the case, for example, in an Fe/GaAs 
contact^ The two-dimensional periodicity of the layered 
systems allows to Fourier-transform the Green function 
in the x and y directions, obtaining a two-dimensional 
Bloch vector k|| = {k x ,k y ) as a good quantum number, 
and retaining an index i to characterize the layer in the 
direction of growth z. The Green function connecting 
the layers i in the left lead and i' in the right lead is then 
written 



VI. THE CONDUCTANCE FORMULA IN KKR 

A. The Green function in the KKR formalism 

In the KKR method the one-electron retarded Green 
function is expanded in terms of local orbitals centered 
at the atomic sites R n as: 

G+(R„ + r,R n , +r';E) 

= -iVEj2RL(r<;E)H2'(r > ;E)S nnl 

L 

+ J2RUr;E)Gli(mL(r';E) (31) 

LL' 

Here, i?£(r;_E) and iJ£ (r;E) are, respectively, the reg- 
ular and irregular solutions of the Schrodinger equation 
for the single potential V n (r) of the nth cell in free space. 
Atomic units are used (e = — \/2, h = 1, m = 1/2). 
The index L — (I, m) represents the angular momentum 
quantum number. The position vector r is confined in 
the atomic cell n; r> and r< are the longer and shorter, 
respectively, of r and r'. In Eq. (|31|l the first term gives 
the on-site contribution to the Green function, while the 
second is the so-called backscattering term, where the 
information on the intersite electron propagation is con- 
tained in the structure constants G^^E). These are 
related to the known structure constants of a reference 
system, e.g. vacuum, via an algebraic Dyson equation 



G+(R + X v + r, R*< + xv +r';E) 
= -J— ( d 2 fc||e* k »(*'^) 

SSBZ JsBZ 

x £ Ri(r; E)G% (k,, ; E)B&, (r' ; E) (32) 

LL' 

where Xv and Xv> are in-plane lattice vectors, R; is the 
interlayer lattice vector, SBZ is the surface Brillouin zone 
of the system and <Ssbz its area. In this equation each 
layer i is assumed to have a unique atom type, hence 
only the index i suffices to characterize the local wave- 
function. In the case of more inequevalent atoms per 
layer, an extra index /i can be introduced to account for 
the propagation between different kinds of atoms. In the 
case of spin magnetism, the Green function is different 
for each spin direction a =\ or J.. The formalism can be 
generalized for fully relativistic calculations, where the 
spin-orbit coupling results in a mixing of the two spin 
channels. 

Once Eq. lj3Tl) or Eq. is substituted into the ex- 
pression l|10l) . with r and r' in different leads, the on-site 
term of the Green function does not contribute, and only 
the intersite term survives; moreover, the spacial deriva- 
tive affects only the local orbital functions R 7 l( r ) an d 
H2 (r), leaving the structure constants untouched. Fi- 
nally, if Eq. (|32(l is used, the conductance appears at first 
k|i-resolved, g(k|| ), which is most convenient in structures 
with two-dimensional periodicity. A k|| -integration then 



7 



gives the result 



Vac 



9 = 



Ss 



BZ JsBZ 



d 2 fc|| 5 (k||) 



(33) 



In Eq. i|l(J|) , both the retarded and the advanced Green 
functions are needed; however, they are related through 
the identity (|18fl for real E. This is used in our formu- 
lation; the energy E is identified in the calculations with 
the Fermi level Ef plus an (in principle infinitesimal) 
imaginary part e which we take very small. 

The case (it) of wire-like structures is completely anal- 
ogous, but the Fourier transform of Eq. Ij32|l is not per- 
formed. More specifically, consider a wire embedded in 
vacuum (see also Fig.|5J|; the vacuum region will be also 
divided artificially in volume- filling cells. One can also 
consider defects within the wire, or even two wires at- 
tached to some cluster of atoms (as indicated in Fig. [21 , 
and solve self-consistently for the electronic structure. 
Due to electronic states at the surface of the wire, the 
first one or two vacuum layers above the wire surface can 
contribute to the conductance, but after that the cross 
section can be truncated. The cross sections left and 
right, where the conductivity tensor must be evaluated, 
consist always of more than one atomic cells, since one 
must include the vacuum region. In this way the Green 
function must be considered for all combinations between 
cells on the left and the right, and the expression for 
the conductance splits up in partial contributions corre- 
sponding to these combinations: 



E E 

/i'(Left) /i(Right) 



(34) 



with g^^i given by Eq. l)10f) but integrated over single 
atomic cells. 

Analogoulsy, in case (Hi) a similar setup can be used if 
we consider semi-infinite 2D leads but current flow local- 
ized in space, as in the case of transport through a con- 
striction. In this case we consider two planes as shown 
in Figure |21 (right), while convergence must be checked 
with respect to the size of the regions considered in the 
summations of Eq. 1341 

In the next subsections we will consider the calculation 
of the spacial derivative of the local radial functions and 
the Green function. 



B. Plane integration 

Firstly we consider a direct evaluation of the conduc- 
tance by use of Eq. I|10|) and calculation of the spacial 
derivative of the Green function at exactly the plane sur- 
faces left (5l) and right (Sr). Both Si, and 5r are as- 
sumed to be in the asymptotic region where the potential 
is stabilized to the bulk one and the evanescent states 
have decayed; in practice one has to perform the calcula- 
tions for several positions of <Sl and <Sr at finite distances 
to verify that the results remain unchanged. 
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FIG. 2: Left: Setup for the calculation of conductance in a 
nanowire configuration. The filled circles represent the wire 
atoms (leads), including defects or an atomic cluster (shaded 
circles). The white circles represent the vacuum region. The 
conductance calculation is truncated at the dashed lines, as- 
suming that the states of the nanowire have decayed beyond 
this region. The arrows labeled z and z' represent the direc- 
tions of growth of the two leads. Right: Similar setup for the 
calculation in the case of nano-size constrictions between in- 
finite leads. Again a truncation at the dashed lines is taken, 
assuming that outside this region there is insignificant tun- 
neling between the leads. 



The set of atomic cells is volume-filling, and the plane 
surface cutting through them inherits the cellular struc- 
ture used in the KKR method; thus it is split in two- 
dimensional tesselating cells. Each one of them is a con- 
vex polygon corresponding to the section of the plane 
that belongs to a convex Voronoi polyhedron^ (or just 
Wigner-Seitz cell in the monoatomic case). In this way, a 
two-dimensional cellular Voronoi construction is defined 
in the plane, each cell of which is completely within some 
three-dimensional cell. Therefore the representation of 
the Green function in terms of local radial functions can 
be readily used. In fact, in systems with two-dimensional 
periodicity, a two-dimensional unit cell consisting of some 
convex Voronoi polygons S^ is constructed and the cal- 
culations can be confined in those. An example of the 
construction for bec (001) surface cells is given in Fig. J3Jl. 

Using such a construction, the k|| -resolved conduc- 
tance is written as: 



s( k n) = -iEE E 

f^fi' LL' L"L"' 

(j£>L»> - Gli (k„ ; E)GIV L ,„ (k„ ; E) (35) 

where we have introduced the KKR current matrix ele- 
ments J^ L ,(m a similar fashion as has been done in the 
paslpi^) in a cell as: 



Jl L ,= [ d 2 r Rl (r; Ep)d z R* L , (r; Ep). 



(36) 



The summation Y] , is over the Voronoi polygons of 
the inequevalent atomic sites of the 2-D unit cells in the 
leads. The calculation of J^ L , is described in Appendix 
IU In the case of a finite system, where a two-dimensional 
Fourier transform is not necessary, the summation is over 
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| S(side) " 

FIG. 3: Two-dimensional geometrical construction for the bcc 
(001) surface cut. The cutting plane goes through the lattice 
sites of an atomic layer (full circles), but must include also 

part of the Wigner-Seitz cell of the next layer (open circles). FIG. 4: Construction for the conversion of volume to surface 
The shaded area shows the two-dimensional unit cell formed, integrals over the lead cross section, 
consisting of two smaller convex Voronoi polygons. 



all atoms in the planes Sl and Sr, and the k|| dependence 
drops. 

The Green function matrix elements G^,, containing 
the information for the propagation from site \j! at a layer 
within the left lead to a site // at a layer in the right 
lead, are off-diagonal, while for the charge density one 
needs only the diagonal (on-site) elements. Nevertheless, 
an efficient O(N) algorithm exists^ for their calculation 
within the screened KKR formalism, i.e., the time needed 
for the calculation scales linearly with the distance be- 
tween the two layers. Thus it is possible to calculate the 
conductance in junctions of more than 100 monolayers 
with present-day computers pS^S 



C. Volume integration: ASA and full cell 

In this subsection we provide an alternative to the 
calculation of the surface-integrated current matrix el- 
ements of Eqs. (|36[1 and (|B4|) . We prove that the cal- 
culation can involve a volume integration over the unit 
cell, instead of a surface one; in principle the results must 
be equivalent, but this method has advantages when one 
wishes to use the atomic sphere approximation. Most 
important is, however, that the /-convergence is much 
better (see Sec. IVIIfl . 

First we observe that the value of the conductance is 
indeed independent of the position of the planes Sl and 
Sr. This can be proven using the fact that the con- 
ductivity tensor is divergenceless (Eq. |7|) . Say that r is 
on Sr; if we consider a second plane surface S R close 
to Sr, we can utilize Gauss's theorem in the volume 
V enclosed by the two planes, to convert a volume in- 
tegral of Eq. in V into a surface integral over Sr, 
S R , plus side-areas. The construction is analogous to 
the one described in Appendix ^ as shown in Fig. 0] 



The contribution from the side-areas vanishes because 
there we have either totally confining boundary condi- 
tions or Born-von Karman boundary conditions leading 
to cancelation from opposite side-areas due to the oppo- 
site surface unit vector orientation; then we are left with 
J s dS<Ty (r,r') ■ z — f s , dS(Tij(v,v') ■ z. The same ap- 
plies for Sl, where z' varies; thus the flux of <Ty (r,r'), 
i.e. the conductance, is independent of the exact position 
of Sl and Sr, q.e.d.. In fact, following these arguments, 
we see that Sl and Sr do not even have to be planes; for 
instance, they can follow the pattern of the Wigner-Seitz 
or Voronoi cells, as long as they meet the requirement 
that they satisfy the Born-von Karman periodic bound- 
ary condition in x and y. In the case that they are not 
planar surfaces, one must of course take the flux of the 
conductivity tensor really along the normal n pointing 
outward at each point of the surfaces, i.e. 

g= I dS I dS'n-a{r,r')-h' (37) 

Since the exact choice of Sr or Sl does not affect the 
result, one can average over the volume V included be- 
tween, say, Sr and S R instead of integrating over Sr. In 
particular, V can be chosen to have a thickness d equiv- 
alent to a unit cell in z-direction, so that one has to 
average over layer-adapted unit cells. In this respect, the 
conductance formula has the same form as Eq. (|35|l , but 
with an extra prefactor of 1/d 2 to account for the volume 
averaging in the two leads (d here is the distance between 
two consecutive lattice monolayers); the current matrix 
elements Jll" , volume-averaged here in the atomic cells, 
have the form of Eq. (|36[) but with the integral being 
three-dimensional over the atomic cell. This can be done 
both in the atomic sphere approximation (ASA) and in 
the full-potential (and full cell) formalism. Details about 
their calculation are given again in Appendix iBl 

A word of caution is due here: it is essential that the 
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ABABA ABABA 




FIG. 5: Two possibilities for cell- averaging of the conductivity 
tensor flux in an ...ABAB... stacking sequence. Left, using 
a cellular division without corrugation of the surfaces, and 
right, using the Wigner-Seitz construction with corrugation. 
The latter is preferable in the KKR method. 

volume averaging leaves no "holes" in the structure. To 
be specific, assume a bcc-like structure. Then, as is evi- 
dent from Fig. [31 the Voronoi cells of the first layer (open 
circles) can touch via the "holes" (corresponding to the 
squares) to the cells of the third layer. In this way, the 
current can partly by-pass one monolayer traveling di- 
rectly to the next one. If one just takes the current 
averaged over the Voronoi cells or the ASA spheres of 
the middle monolayer, one forgets to calculate this part 
of the current; this is why the full many-atom unit cell 
must be taken, so that no such holes are left. 

If one uses the usual Wigner-Seitz cells (or atomic 
spheres), Eq. I|B5|) can induce a small inaccuracy. The 
reason is that the volume constructed by such cells is 
in general not included between planar surfaces, as in 
Fig. ijlfl. but rather between corrugated surfaces, in ac- 
cordance to the form of the Wigner-Seitz cells. In such 
a case, the full conductivity tensor er(r,r') (not just the 
cr z z component) and Eq. I|37|l should be used in princi- 
ple. To avoid such a more complicated calculation, two 
ways can be followed, as demonstrated in Fig. [5] Firstly, 
one can persist in using layer-adapted unit cells (par- 
ellelepipeds), which give no corrugation. This has the 
disadvantage that such cells can be too flat so that the I 
expansion of the cell-centered KKR Green function and 
wavefunction converges poorly. Secondly, one can aver- 
age over more than one monolayers; in this example they 
would be A-A, A-B, B-A, and B-B. Then the corrugated 
region is a smaller fraction of the total averaging volume, 
so that the error becomes smaller. Test calculations on 
this will be given in the Section IVlII for bulk Al. 

D. Current matrix elements and selection rules 

Which method, volume or plane integration, is most 
convenient and accurate depends on each specific prob- 
lem; however one can have "rules of thumb" on the diffi- 



culty and convergence of each one. In many systems the 
atomic sphere approximation is used, where the poten- 
tial around each site is assumed to be spherically sym- 
metric, but still a full multipole expansion of the charge 
density is taken. This has the advantage of greater sim- 
plicity and less computational effort than a full poten- 
tial and full cell description. In such cases, the plane 
integration is not applicable, since a plane cut through 
volume-conserving spheres cannot give an accurate sur- 
face area; the volume integration is then the only way. 
In the case of a full cell treatment, when one has the 
correct Wigner-Seitz or Voronoi volume tesselation, one 
must consider that the plane method has a drawback, 
namely that the plane might go through regions only at 
the edge of certain cells, where the Z max cutoff seriously 
affects the accuracy of the results; on the other hand the 
volume integration averages out such inaccuracies. 

One has also "selection rules" that make certain J 
elements vanish. This is most easily seen if one uses 
spherical potentials. To be specific, say that the plane 
goes through the atomic site at z = 0; then, in the plane 
integration one can easily see that the elements Jll> are 
non- vanishing for I' = i, 3, 5, . . . when / is even and for 

I' = 2,4,6, when I is odd. On the other hand, in 

the case of volume integration the Jll' are non- vanishing 
only for I' = l,l± 1, i.e. Jll' is band-diagonal in I and 
/'. This can be viewed as an advantage of the volume- 
averaging method, since it means that, if one describes 
the electronic structure with orbitals truncated at Z max , 
for the accuracy of Jll> one has to consider wavefunc- 
tions only up to I' = ^ max + i. A full cell treatment adds 
more nonzero elements, but the main contribution still 
comes from the ones mentioned. 



VII. EXAMPLES 
A. Band counting in bulk conductance 

When the Landauer formula is applied to a perfectly 
periodic material, e.g. the bulk of a crystal, it gives a 
finite conductance which physically represents the con- 
ductance of a long wire placed between two phase- 
randomizing electrodes^ Resolved in k||, the value of 
the conductance equals the number of right-propagating 
(or equivalently left-propagating) states at Ep for this 
k||. In other words, one has to count the Fermi sur- 
face bands for that k|| , which propagate in the direction 
k z _L ky with v z > 0. This is demonstrated in Fig. |SJ 
where part of the Fermi surface of Al is presented, in 
the k x -k z plane, together with the conductance in the z 
direction as a function of k x (with k y — 0). Actually, 
some of the bands shown have v z < 0, but their equiva- 
lents with v z > exist symmetrically for k z < 0. Clearly 
the conductance (in units of e 2 /h) equals the number of 
bands at Ep for each k x , giving a stepwise picture. 

Also in Fig. we can compare the results for angu- 
lar momentum truncation at Z max = 2, 3, and 4. In- 
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creasing Z max results in a more accurate description of 
the wavefunction derivative. Basically the slower conver- 
gence here arises from the relatively high Fermi energy of 
Al, and is also present in a free-electron approach for high 
Ep. As noted in the previous section, if the wavefunc- 
tion is accurate for some l ma .x, for the derivative one has 
to take Z max + 1. For most practical purposes, / max = 3 
is enough, considering also that the calculation time for 



4 is more than three times the one of L 



3 (due 



to matrix inversion, the calculation time scales roughly 
as (Z max + 1) 6 )- 

In Fig.Qthe conductance for the same system is shown, 
but also analyzed in the various interlayer contributions. 
We see that these can exhibit fluctuations in k|| , reflecting 
fluctuations as a function of the interlayer distance. How- 
ever, when averaged over two monolayers on the left and 
two on the right, the fluctuations practically cancel each 
other. The origin of these fluctuations is the Wigner-Seitz 
construction for the unit cell, resulting in corrugation of 
the surfaces where the current is calculated. Due to the 
fact that we account only for the zz-component of the 
conductivity tensor, when we have corrugation the ma- 
trix elements in Eq. (|15fl arc not integrated correctly and 
the nondiagonal current matrix element do not vanish; 
thus, beating effects of the conductance appear. As the 
conductance is averaged over more than one monolayers, 
the corrugation-free region in the middle increases, the 
relative error due to the corrugation decreases and the 
steps in the conductance become flat. If one would use 
tetragonal or prismatic unit cells, suited for layered ge- 
ometry, then the fluctuations would be absent; however 
such cells are not well-suited for an accurate calculation 
of the wavefunctions in KKR. 

Finally, in Fig. [S] we see a calculation for the same 
system, but employing the in-plane integration, rather 
than the volume averaging. Here the convergence with 
( ma x is poor for the reasons explained in the previous 
section; even for Z max = 4 the deviations from integer 
conductance values are large. 




0.2 0.4 0.6 0.8 
kr 



FIG. 6: Fermi surface in the k x -k z plane (bottom) and con- 
ductance (in units of e 2 /h) as a function of k x (top) for bulk 
Al and for (Z max = 2, 3, and 4). 



B. The effect of the nonzero imaginary part of the 
energy 

Although the conductance should be calculated at a 
real energy E, the Green function in the KKR method is 
always calculated at a complex energy E+ie, and the real 
E is approximated by taking very small, but nonzero, e. 
This can have an artificial damping effect to the conduc- 
tance, since waves within a small energy range around E 
are effectively superimposed and finally the phase is ran- 
domized, especially if the leads are seperated by a large 
distance. In Fig. [5] we show an example of how small 
e should be in a realistic calculation. The system here 
consists of two Fe leads with parallel magnetic moment 
seperated by a ZnSe spacer. Electrons are injected from 
the first lead into the ZnSe conduction band, and are 
detected by the second Fe lead. The spacer thickness is 



varied from 9 to 97 monolayers, and due to the multiple 
reflections at the two interfaces transmission resonances 
appear at certain thicknesses. For more information we 
refer to Ref. |3(]. The conductance of the majority elec- 
trons for k|| = 0, i.e., at the f-point, is presented in 
Fig-EDfor a choice of e = 0.02 mRy, 0.2 mRy, and 1 mRy. 
For e = 1 mRy there is strong artificial damping, while 
e = 0.02 mRy is adequate even for the large thickness 
of 97 monolayers. Note that this damping cannot model 
an effect of temperature, because, when we depart from 
the real axis, the spectral density of a state transforms 
from a delta function to a Lorenzian distribution and not 
to the derivative of the Fermi function. Due to the long 
tails of the Lorenzian the damping is much stronger than 
for a Fermi distribution of the same halfwidth. 
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FIG. 7: Contributions to the fee Al bulk conductance from 
the various layers (A-A, B-B, B-A, and A-B in the ...ABAB... 
stacking sequence of Fig. |K| right) (dashed lines) and aver- 
aged sum (full line), as a function of k x , for / max = 4. The 
fluctuations in the partial contributions come from artificial 
beating effects due to corrugation, and vanish in the averag- 
ing procedure. The inset shows some of the fluctuations in 
more detail. 

VIII. SUMMARY AND CONCLUSIONS 

We have presented a formalism for the calculation of 
ballistic conductance in solids, based on the KKR Green 
function method for the ground-state electronic struc- 
ture combined with the Landauer-Biittiker approach. It 
makes use of the result of Baranger and Stoned con- 
necting the derivative of the one-electron Green function 
to the conductance. For the foundation of the formal- 
ism, we have discussed the relation of the S'-matrix be- 
tween Bloch in- and out-states to the conductance. We 
have given an expression connecting the S'-matrix to the 
Green function of the system, generalizing the theory of 
Baranger and Stone to include the realistic band struc- 
ture of the leads. 

The convergence of the method with angular momen- 
tum cutoff (/ m ax) was studied and found to be compara- 
ble to that of KKR. It can be applied to systems with 
two-dimensional periodicity as well as nanowires. Our re- 
sults show that the volume integration and averaging of 
the current matrix elements, applicable both in ASA and 
full-cell or full-potential approaches, gives well-converged 
results for the calculation of ballistic transport. Owing 
to the linear scaling of the calculational effort with the 
number of layers of the screened KKR formalism (O(N) 
scaling), our method is suitable for large systems. 



FIG. 8: Conductance as a function of k x (top) for bulk Al 
(/max = 2, 3, and 4) in the case of plane-integration. It is not 
necessary to average over more monolayers, but the I conver- 
gence is rather poor and the integer valuer are not reached 
even for l m&K = 4. 
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FIG. 9: The effect of the nonzero imaginary part e of the en- 
ergy. Conductance (in units of e 2 /h) as a function of spacer 
thickness for spin injection through the conduction band in a 
Fe/ZnSe/Fe (001) junction. The oscillations are due to multi- 
ple reflections. For e — 1 mRy there is strong artificial damp- 
ing, while e = 0.02 mRy is adequate even for the thickness of 
97 monolayers. 



Acknowledgments 

The authors are grateful to Professor N. Stefanou for 
helpful and motivating discussions. Moreover, financial 
support from the RT Network of Computational Mag- 
netoelectronics (contract RTN1-1999-00145) of the Euro- 
pean Commission is gratefully acknowledged. 



12 



APPENDIX A: CURRENT MATRIX ORTHOGONALITY RELATIONS 



Let and ^k'a' be two Bloch wavefunctions of the same hamiltonian at the same energy. Then we shall prove 
that the following relation holds: 

J dS(W ka; v a ,) z := J dS(*tJv z Vv a >) = ^|(«ka)*|<W<W, (Al) 

where S is an (infinite) cross-sectional area in x and y directions. 

The proof has as follows: First we note that, as a consequence of the single-particle Schrodinger equation for a real 
potential, 

^ — ^ 2 7 77, 

VW ko;kv = V(* kQ V tf k , a ,) = - — {E k , a , - E ka )*£ a * k , a , = 0, for E k , a , - E ka (A2) 

where the band indices a and b are used explicitly. This is just an expression for current conservation of Hamiltonian 
eigenstates. Then, for each volume V enclosed by a geometrical surface S, Gauss's theorem gives 

dSh • W ko;kV = / AVW ka;kV = (A3) 
s Jv 

where h is a unit vector at the surface S pointing outward. In particular, V can be chosen as a prismatic normal 
cross section of the lead, extending from z to z + d. Then its surface S can be decomposed in two plane cross sections 
Si at z and S2 at z + d, as the bases of the prism, plus side-areas S'side at the lead surface, as shown in Fig. 0] At 
these side areas we either have confining boundary conditions, i.e. ^kalside = 0> whence W ka;k / a / | s id c = or Born - 
von Karman periodic boundary conditions, whence for each prism side there is the opposite one with the same value 
of ^kalside and W ka;k ' / | s id e but with opposite orientation unit vector n; then the sum of their contributions to the 
surface integral will be again zero. In this way, we are left with the two bases of the prism; they have opposite unit 
vector orientations, thus: 



/ dSi(-z) ■ W ka;k / Q / + / dS 2 z ■ W ka;k / a / = 

JSx J S 2 

J dS , i(W ka;k / a /) z = / dS , 2(W ka;k / a ') z 



(A4) 



(A5) 



i.e. J s c?5'(W ka;k / a ') z is independent of the position z of the cross-section S. This means that one can average it over 
the whole Born - von Karman supercell (in z) of length L and volume VbvK = S x L: 

[ dS(W ka;k , a >) z = ~ I dz ( dS(Wka;k'a')z(x,y,z) = y { d a r(W ka:k , a ,) z (x,y : z) (A6) 
Js L Jo J L 7v BvK 

We now convert the volume integral in a sum over unit-cell integrals J v d 3 r and employ the Bloch property 

W ka;kv (R, + r) = e j(k '- k)R *W ka ;k'a'(r) (A7) 

to get: 

/ dS(W ka;k ,Mx,y,z) = )E/ d 3 r(W ka:k -^) 2 (R ! + r) (A8) 

JS L • JVg 

N f 

= -pSkk' / d 3 r(W ka:k 'aO z (r) (A9) 

L JVa 

where N is the total number of lattice sites in Vb v k, and 

e i(k'-k)R.< = ^ kk , (A10) 
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has been used; Vq is the unit-cell volume. Note that NVq = VbvK = SL, thus N/L = S/Vo, and we get: 

[ dS(yf 1at;Va ,),(x,y t z) = £r6 kk , [ d 3 r(W ka:kW ),(r) (All) 

JS V JVn 



For k k' this gives zero, while for k = k' and a = b we have 

/ d 3 r(W kaM Uv) = ^ / d 3 r(*L — V* ka - * ka — V^J = i^v ka (A12) 
Jvo n Jvo im im h 

Eqs. I|A12() and IjAllI) verify (|A1|) except in the case of band crossing, when k = k' but a ^ a'. In this case we may 
use the identity connecting the momentum operator p op to the Hamiltonian H and the position operator r op : 

ITYI 

Pop = -^[H,r op ] (A13) 



For the evaluation of (| A6|> we need the matrix element of p op := (h/i)W. Using (|A13|I we get 

h 

'^BvK ^ VfivK 



( d 3 r^i aPop ^ k , a , = (E ka - £W) f d 3 r^ Q r* k ^ (A14) 



In case of band crossing, (i? ka — E^i a i) = 0, but this does not mean that the whole expression vanishes, since the 
integral might diverge. It can be calculated in a standard way by utilizing Bloch's theorem and reducing it to the 
unit cell. We have: 

/ d 3 r** Jr)r<I/ k , a ,(r) = W dV**^ + r) (R. t + r) ¥ k , a ,(Ri + r). (A15) 

JV BvK i JVo 

Using the Bloch properties of ^ ka and 'JVq' and the relations: 

WCtf-ion, = (ML S (k'-k) (A16) 

J2 e i(k'-k)R iR . = V k /e i ^- k ) R ' = ^V k ^(k' - k) (A17) 

i i 

we obtain after some manipulations 
A f rf 3 r* ka Po P * k ' Q ' 



^-5(k' - k) !(E ka - Ex*) I d 3 r^* ka r^ k , a , ~ hw ka S aa , - (E ka - £W)V k [ d 3 r* kQ * k , a , 1 
Vo L Jvo JVo ) 



(A18) 



In all terms, k = k' can be directly substituted due to the ^-function, except in the last one, where the one must first 
perform the integration and the derivation. In the second term, the orthogonality relation J v d r\I/ k / a \I/ k v = ^aa' 
after the substitution k = k' has been used. From this expression we immediately see that in the case of band crossing, 
i.e. k = k', Eka = S k v> but a ^ a', the expression vanishes, so the proof is complete. In passing we note that, if 
a = a' , the expression gives the group velocity as expected. 

APPENDIX B: CURRENT MATRIX ELEMENTS IN KKR 

In the plane-integration formalism the KKR current matrix elements read 

Av=\ d 2 rR L {T-E F )d z R* L ,{v-E F ) (Bl) 

J Si 

where Si is the surface cut of the atomic cell i with the plane passing through it. In the full-potential KKR formalism, 
the wavefunctions are expanded in terms of real spherical harmonics as 

RL(r) = J2l R ^Ur)Y Ll (9,<t>) (B2) 

Li r 
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The real spherical harmonics are of the form 

Y L (0, 

where 

CtL 



a L ■ P;' m '(cos0) • trg(m0), 



(B3) 



121 + 1 (l-\m\)\ 
2tt (I + \m\)\ 



22 + 1 

47T 



, m = 



P ; ' m '(cos#) are the Legendre functions, and trg(m</>) = costo.0, if to > or sin | , if to < 0. Thus one has to 
decompose d z into d r , dg and d<f,. The first affects only the radial part Rl 1 l(t) and is calculated numerically; the 
other two affect only (0, <j>) and are calculated analytically. After some algebra one arrives at the result 



Jll' — 



dr — 



J> ia [(rd r Rl 2L ,(r) - (l 2 + l)Rl 2L ,{r)) uP 1 ^ 1 (cos 0) + (l 2 + KDPfcWj^W 

r<i> ] t 

2_^a Ll P}" lll (cos9)R LlL (r)2_^ / trg(mi^)trg(m 2 ^) # . 



(B4) 



Here, r m ; n and r max are the radii of the inscribed and circumscribed circles, respectively, of the convex polygon, on 
which the (^-integration is performed, with center the z-projection of the atomic site on the plane; Gentry an d '/'exit are 
respectively angles of entry into and exit from the convex polygon as the ^-integration is performed. 
In the volume-averaging formalism, the KKR current matrix elements have the form: 



Jlv = / d 3 rRl(v)d z R^(v) 

J coll 



(B5) 



These can be computed within the full cell or ASA formalism; here we shall present both results. 

In the full cell approach, the potential is truncated at the boundary of the Voronoi atomic cell. This is achieved by 
introducing the characteristic, or "shape", functions 9(r), being equal to unity in the cell and vanishing outside4£ 
Their expansion in spherical harmonics, 



e(r) = ^e L (r)y L (#,, 



(B6) 



is used in the calculation of the current matrix elements. After some manipulations we obtain 



Jll' — 



d\Q{v)R L {v)d z R* L ,(r) 



(B7) 



\vs 



EEE [ / dr (RL lL (r)d r Rl 2L ,(r) ~ l l±lR LiL ( r )Rl 2L ,(r)) G L3 (r) J- £ C LxL2Ll C x ^ hi 

r.- r._ r_ L/ \ T / r . 



h\ hi L 

i 

r 



dr ^R LlL {r)R* L2L ,(r)Q L3 (r)-^^—(l 2 + \m 2 \)Ct lLa h-i,m3 

Oilo — l.mo 



(B8) 



where Cl 1 l 2 l 3 = J dflYi ll (f)Yi l2 (f)YL 3 (f) are the Gaunt coefficients and the identity J dQY^Y^Y^Y^ = 
C LlL2Lll C LoL3Li has been used. 
The ASA result is simpler since it does not involve the shape functions. The local orbitals have only a spherical 
part, 



r l (v) = -ij,(r)y L (n) 

r 



(B9) 



whence the current matrix elements become 

1 

"1,0 

, {l' + \m'\)ai 



Jll< 



Cluxo / drRi(r)d r Rf,(r) 



I' + l 



ai.o 



Cll>xo / dr-Ri(r)Rp(r) 



Oil'-l, 



-l<5n 



dr-Ri{r)R* v {r) 
r 



(BIO) 
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It is implied that the integrals are within the atomic sphere. 
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